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Abstract 



0^ , This paper extends the work of Clarke |l| on the Bayesian foundations of the biomagnetic inverse 

5»j ' problem. It derives expressions for the expectation and variance of the a posteriori source current prob- 

C/3 ' ability distribution given a prior source current probability distribution, a source space weight function 

y , and a data set. The calculation of the variance enables the construction of a Bayesian test for the ap- 

C/3 ' propriateness of any source model that is chosen as the a priori infomation. The test is illustrated using 

^~>| both simulated (multi-dipole) data and the results of a study of early latency processing of images of 

i-Q . human faces. 
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1 Introduction 

The magnetoencephalographic (MEG) inverse problem (sometimes known as the biomagnetic inverse prob- 
lem) is the process of identifying the source current distribution inside the brain that gives rise to a given 
set of magnetic field measurements outside the head. The problem is difficult because the detectors are 
limited in number and are sensitive to widespread source currents, and because of the existence of silent 
and near-silent sources. 'Silent sources' are configurations of current density inside the brain which give zero 
magnetic field outside the head (e.g. all radial sources are silent when the head is modelled as a conducting 
sphere). It follows that the general biomagnetic inverse problem is ill-posed and under-determined. 

The most common way of reducing the problem and rendering it tractable is by characterising the 
source in terms of a limited number of effective current dipoles. Such source descriptions provide links 
with the dominant functional architecture model of the brain in which processing is described in terms of 
localised activity with interactions between essentially separate regions. Multiple dipolc models have enjoyed 
considerable success in the analysis of sensory and motor cortex (e.g. |g|, g, ^, ^). 

Growing evidence for the existence of more diffuse brain networks have led to an interest in distributed 
source algorithms. Several have been proposed [0, M, R H nft O. These algorithms have been designed to 
cope with the non-uniqueness of the problem, primarily by restriction of the source space and by regularisa- 
tion. Each algorithm leads to a unique solution (from the infinite number available) through its particular 
choice of source basis, weight functions, noise assumptions, and, in many cases, cost function. There has 
been an extended debate about the accuracy and value of these methods. This proceeds at two levels; the 
technical ability of the various algorithms to recover a simulated source distribution (often quoted in terms 
of one or more source current dipoles), and the electrophysiological appropriateness of the type of source 
structure favoured by particular algorithm parameters. So, for example the simple minimum norm solution 
m , which tends to produce a grossly smeared and superficial source distribution may be compared with the 
LORETA solution fTq] which favours smooth but regionally confined current distributions. The issues have 
been fully debated in recent conferences |]1^ [l3|). 

The many to one nature of the mapping of sources to magnetic fields suggests that a probabilistic approach 
to reconstructing the sources from the magnetic field could be used. A Bayesian probabilistic approach was 
first proposed by Clarke |l[ . More recently, Baillet et al have described an iterative approach which combines 
both spatial and temporal constraints within a unified Bayesian framework designed to allow the estimation 
of source activity that varies rapidly with position, e.g. across a sulcus [|lj]. Schmidt et al have developed 
a probabilistic algorithm in which a bridge is made between distributed and local source solutions through 
the use of a regional descriptor within the source representation ||l^. In this case, a Monte Carlo method is 
used in the absence of an analytic solution for the expectation value of the source current. 

Here, we are proposing an alternative Bayesian approach. It includes the explicit assumption of both a 
prior source current probability distribution and a source space weight function, and allows the calculation 
of the expectation and variance of the a-posteriori source current probability distribution. The derivation 
of these quantities is detailed in Section 3. The inclusion of the prior probability and the calculation of the 
variance provide a means by which the consistency of a model (assumed as the prior) can be tested against 
the data to reveal those parts of the source distribution that are statistically robust and, conversely, where 
the model is inadequate. Numerical calculation of significance is possible. A straightforward extension of 
this idea is the direct comparison of two data sets to reveal where, within a given model, there are significant 
differences in their associated source distributions. In the final part of this paper, both simulated and real 
data will be used to illustrate these various uses of the technique. 

2 Specification of the problem 

The arrangement of sources and detectors for the biomagnetic inverse problem is shown in Figure |l| The 
sources giving rise to the measurements are assumed to be restricted within a source space Q, which may 
be smaller than the whole head volume (e.g. if the sources are assumed to be cortical). The current density 
within the volume Q is assumed to belong to the space of square integrable vector fields on Q, which we call 
J. 

The measurement process typically gives successive sets of data (~ 100 channels) every millisecond. In 
this paper the data for each time instant is processed independently, and the data from a single time instant 




Figure 1: A schematic experimental geometry. 

is collected into a vector m G M . If j(r) € J then the measurement process can be represented by a 
functional z : J -^ R .A subscript notation will be used to identify the sensor, i.e. Zi is the ideal reading 
from the ith sensor. So the basic equation is: 

rht = Zi{j) + Ci (1) 

where the e^ are the measurement errors, which are assumed to be normally distributed with zero mean 
and covariance matrix a^D (where a is the standard deviation of the errors and D is a symmetric, positive 
definite matrix). 

To compute the functional z{ ) on a computer (i.e. to solve the forward problem) requires a volume 
conductor model of the head. In this paper the precise model used is irrelevant, so the final results will be 
written in terms of z. This is done is via the Green's functions Li for the problem, which are defined by 



z,(j)- / L,{r)-jir)df (2) 

JQ 

Stated simply, the inverse problem is to estimate j{f) given the data vector m. Obviously the given data 
are not enough to determine j{f) uniquely (for several reasons). The approach adopted here starts from the 
same point as used in Clarke 0|, a statement of Bayes's theorem: 

POeA|„.ei,,.S2il|lL^_5ff(il^ (3) 

V(m E B) 

where A is a set of currents and B is a set of measurements. This equation reads, the a posteriori probability 
of a current set A after the measurement B is proportional to the probability of producing the measurement 
B given that the current is in the set A times the a priori probability of the current set A. {'P{m £ B) is a 
constant for any measurement set B). 

In this paper the probability distributions (both the prior probability and the errors) are assumed to 
be gaussian and so it is permissible to work with probability density functions. A further simplification 
is achieved by shrinking the measurement set _B to a single point {m} (this ignores the finiteness of the 
precision of the measurements). Equation ^ then becomes: 

PmO) fx p(]) X e(m - z(J)) (4) 

where p is the a priori distribution, p^ is the a posteriori distribution and e is the error distribution. 
Throughout the paper, probability density functions will only be determined up to a constant. The constant 
of proportionality is found by requiring that the probability is normalised to one. In this paper both p and e 
are assumed to be Gaussian and then p^ is calculated to be Gaussian. An error probability density function 
e consistent with the Gaussian assumption may be written: 

e(e)aexp|-2^e-i?-ie| (5) 

This generally valid expression will be retained throughout the derivation in this paper. In practice, the 
noise covariance matrix D may be difficult to estimate and, for simplicity, the simple form D = I will be 
used in the later illustrations. 



3 Derivation 

In this section formulae for the maximum hkehhood current distribution and also the expected error distri- 
bution are derived under specific assumptions. First, an inner product on J is defined by: 

where cij(r) is a weighting distribution defined on the source space Q. This provides a method of inputting 
prior information of the location of sources (e.g. gained from MRI images) into the algorithm. 

Clarke |l) assumed that the maximum likelihood prior current density was identically zero. Here that 
restriction is avoided and an arbitrary prior current j^ will be introduced as a parameter of the method. 
The a priori probability distribution on J is defined using j^ and the inner product: 

p(/) (X exp |-^ (J- JP , J " J^)} (7) 

where (3 is the assumed a priori standard deviation. 

To proceed further a basis is needed for J. A 'natural' choice is a basis that is related to the measurement 
functional z. So an obvious candidate is a basis derived from the adjoint map to the measurement map from 
J to R . This gives a map z : M — > J (since J is self-dual) defined by: 

z{a) , j) = a^zij) (8) 

Explicitly this gives the set of linearly independent distributions {Lu{f)Li{f)}. This set is extended into 
a basis of J that includes the silent sources by adding vectors {Ci} which are chosen to be orthogonal to the 
{luLi} i.e. 

(c^i. ,4)=0 Vz,j (9) 

Since {u!{f)Li} [J{Ci} is a basis of J a general current density j(r) can be written in terms of this basis as 

N 

i=l i 

To simplify the notation the components of currents are written in column vector notation: 

A simple computation gives 



where Py = (ujLi , uiLj) and Qij = (Ci , Cj) since by construction (ojLi , Cj) ~0. 

Now the two a priori probability density functions (Equation M and Equation H) may be combined with 
Bayes's theorem (Equation H) to obtain the a posteriori probability density: 



a-aP\ fP 0\ fa-aP\\ f 1 ,. ,^..^^-i.~ 



p^(i)(xexp|-^(^^_^pj (^0 gj I^6-6PJ/''''P\"2^(™"^(-^')) -^ (^-^(j))| (13) 

The task now is to manipulate this equation so that it is explicitly in the form of a gaussian distribution. 
As a first step the exponentials are combined: 



p^(y)ocexp -— ^ 



1 

2^2 
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(14) 



where C = c? j (? and z{i) has been replaced by Pa. Next, the terms involving operators on a are simplified by 
completing the square All constant terms can be absorbed into the normalization constant of the probability 
density function and are ignored in this derivation. 

Operators on a = (a - aFf C,P{a - aF) + (m - PafD-^{m - Pa) (15) 

= a'^CPa - 2a'^CPaP + a'^PD-^Pa - 2a'^PD-^m + const. (16) 

= a^ (CP + PD-^P) a - 2a''[PD-^m + (.Pa^] + const. (17) 

= {a- hf (CP + PD-^P) {a-a) + const. (18) 

where a is defined by: 

(CP + PD-^P) a = PD^'^m + (Pa^ (19) 

So, a ^ (CD + Py^rh + QDa^) (20) 

A modified expression for p^{j) is now available, 

p^(J)ocexp|-^r:Jj K^ ^^^^ Iq]\1_Ip]\ (21) 




I I (a-a\ ((PD-^P + CP) " \ (" "\i ^oo^ 

°'"^P|"2^(,6-6^j (,0 — 1-) (.-.pK (22) 

This is explicitly in the form of a gaussian from which the mean current and covariance matrix can be 
identified by inspection. At this stage it is clear that the mean value of b is V i.e. that there is no change 
from our prior knowledge (this is because by construction all the information that the experiment provides 
is orthogonal to the Ci). 

In order to produce images of the a posteriori current density, it is necessary to find the distribution of 
a single statistic that can be computed. The statistic is defined through a 'test current' t — xv^i^f^a where 
XVk {^ is the characteristic function of a voxel in the brain and Ca is a unit vector. This is another departure 
from Clarke [y in which a delta function test current is assumed. This choice causes problems because the 
inner product (t , t) which is needed below (see Equation Bq) is undefined for a delta function. 

The distribution of the statistic A = It , j) will now be determined. First, the coefficients of the basis 

elements required to construct t are identified: 



A = (i , j) (23) 

= /r, ^a,c^(r)i,(r) + ^6.A(r)\ (24) 

= ^ a, (r, uir^Uf)) +J2b^ (r, A(f)) (25) 

i i 

= u^a + v^b, say. (26) 

Equation G3 is projected onto this particular linear combination of co-ordinates to find the probability 
density of A 

,., f 1 {u-'a + v-^b-u'^-a-v-^bPf \ _. 

p™(A cxexp<^ -— ^ ^ . '- ) (27) 

™ \ ^a^u-^{PD-^P + CPy\ + v-^Q-^v/(:\ ' 

The mean of A can be identified from Equation ^ by inspection. 

mean of A = u'a + v'^b^ (28) 



The term v'^b^ cannot be computed directly since the basis fields Ci have not been defined explicitly. The 
problem may be overcome by expanding (t, j^ 



*^' ^7 = y'T. <^i^U^ + E ^'^^^^) ) (29) 

= ^ < (r, u^{r^u^) + Y,bi (r, A(f)) (30) 

■i i 

= u'aP + v'^bP (31) 
Equation ^ may now be rewritten using only references to known vector fields as: 

mean of A = u'h + (f, ]p) - u^a^ (32) 



u'{C,D + Py^ (m + CDaP) + /f, f\ - u^a^ (33) 

^^((1? + Py^ {m - PaP) + (t, jP) (34) 



This equation explicitly writes the expectation value of the statistic X — (t , j) as a sum of two terms. 
The second term is the statistic for the prior current, and so the first term can be identified as the correction 
to the prior suggested by the measurements i.e. the first term shows the difference between the expectation 

oi (t , j) before and after the experiment was made. This is the first central result of this paper and it is 
worth stating explicitly: 

Change in expectation of (t , ) = u'^{C,D + Py^ (m - z{f)\ (35) 



Using Equation [2^ it is possible go further than this and determine the statistical significance of the 
statistic. This is because the variance of the variable A can also be read off from Equation p^ as: 



variance of A = a 



M^(PL'-ip + CP) ^u+jv'^Q'^v 



(36) 



In order to derive an expression in the form of computable matrices, the term v^Q ^v must be rewritten. 
To do so, t is written as a linear combination of basis elements: 

t{r) = J2 ^.^ir)U (r) + ^ 2/< A (r) (37) 

i i 

Of course x and y are related to u and v. In fact: 

u, = (r, u;{r)Ur^) (38) 



= (y2 ^M^L, (f) + ^ y,C, (f) , ^(f)i, (f) \ (39) 

= ^ X, (cj{f)L, (f ) , u; {r)L, (f) ) (40) 

j 

= Y.'^JP^, (41) 

3 

Similarly, v = Qy. Using these relationships, the inner product of t with itself can be computed: 

(t' ^ = r£xM^U{r^ + Y. y.C,{r) , ^ x,a;(f)4(f) + ^ 2/,4(f)\ (42) 

\ i i j j I 

= x-^Px + y^Qy (43) 

(44) 



from which, 

(t, i) =u'^p-^u + v''Q-\ (45) 



Equations |36| and ^ can now be combined to generate the following formula for the variance: 

(46) 



variance of A 



i'{PD-^P + C,P) ^u+j{{t,T)-u'p-^u) 



This equation is a generalisation of the results of Clarke M , some consequences of it were explored in [nq 
and p7| . It is interesting to note that the existence of a prior current density does not affect this variance. 

It may be helpful to relate features of this equation to the measurement system. The second term in 
Equation Efl is multiplied by a^/C = /3^ and so is independent of the assumed noise levels in the detectors. It 
represents a variance derived from the finite number of measurements and the geometry of the experiment. 
Since it is proportional to [3'^, it can be interpreted in terms of the truncation error becoming less and less 
important as the certainty of the prior distribution jP increases. 



The first term in Equation 46 is proportional to a^. It shows how the noise in the data is reflected into 
source space. The unregularised form of this term was derived previously by loannides et al in |18| using 
an ad hoc argument. loannides et al obtained the regularised form by replacing occurrences of P~^ in the 
unregularised form by {P + ^/) . In the notation used here this gives a first term as follows 

first term = w^(P + Ciy\P + (iV^u (47) 

This does not agree with Equation ^ which has the form (in the case of the measurement error being 
uncorrelated, i.e. D = I) 

first term = u^p-^{P + C,iy^u (48) 

The behaviour as C tends to infinity is that the variance tends to zero. This is reasonable since, for fixed 
experimental noise levels (i.e. fixed a), C, tending to infinity corresponds to /3 tending to zero, which in turn 
corresponds to greater and greater certainty that the prior is correct. When j3 is zero the a priori current 
distribution is known with absolute certainty. This is consistent with the above analysis which indicates 
that, in this case, the a posteriori current density is certain to be equal to the a priori current density. 

Note that, when computing the term v7 P^^u in Equation E6^, any reasonable algorithm (e.g. Choleski's 
algorithm) for computing P^^u can be used even though the matrix P is ill-conditioned. This is because 
the large residual vector which results in this calculation is annihilated by the inner product with u. So the 
computation of the whole term is well conditioned. 

4 Applications 

The main analytical results of this paper (Equations ^, |3^ and ^) provide the means of solving the MEG 
inverse problem with specific assumptions and of assessing the robustness of the solution. In this section, 
this approach will be illustrated through three studies: a simulation of a few-dipole source set; an analysis 
of the appropriate dipole model for data from a real experiment on face processing; and a comparison of 
responses to different visual stimuli from the same real experiment. 

All the illustrations are based on the same experimental arrangement and the same instrument, the 
Neuromag-122-^^'^ ]1^. This is a helmet MEG system that contains 61 pairs of first order gradiometers 
{dBr/dO, dBr/d4> in spherical polar co-ordinates), covering the head (Figure ^. The outputs of each pair 
of detectors are closely related to the dominant tangential ionic current fiow in the region underlying the 
relevant sensors. Also shown in Figure g is the assumed source space, a 2-d spherical shell of radius 0.08 m 
covering a 2 radian by 2 radian solid angle over posterior regions of the cortex. 

The source configuration for the first simulation study was three dipolcs distributed in an isosceles 
triangle configuration (co-ordinates (0,-0.08 m,-0.02m), (-0.05m,-0.06m,0.02m) and (0.05m,-0.06m,0.02m) 
with orientations (1,0,1), (1,0,-1), (1,0,1) respectively). The positions of the dipoles were chosen so as to 
represent approximately a central primary and two bilateral secondary visual processing areas. The precise 
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Figure 2: The experimental geometry, a) and c) show the 61 measurement sites in the helmet arrangement. 
Each square represents a pair of orthogonal detectors. Also shown is the source space, consisting of a 
discretised mesh covering a part of a spherical shell, b) Shows the coordinate system with the x-axis along 
the line joining the pre-auricular points and the y-axis joining the inion and nasion. 



locations are not exactly on the source space shell but are displaced by 3 mm, 1 mm, and 1 mm respectively 
from the shell. The dipolc locations and activation curves are shown in Figure 0. The central source is 
activated first, followed by synchronous activation of the lateral sources. The forward problem is solved 
using a homogeneous sphere conductor model centred at the origin. Gaussian noise has been added to the 
computed dipole signal so that the integrated noise power is equal to 50% of the integrated signal power. 
Examples of simulated data are shown in Figure 0c. 

The simplest approach to the inverse problem is to use Equation p4 with a zero prior current distribution. 
This simplification results in the same formulae as the probabilistic algorithm that has been used for several 
years pO| , pl[ . The resulting expectation of the a-posteriori current distribution is shown in Figure 0a. 

However, using our analysis, it is possible to employ an approach which goes further in comparing 
different source descriptions. Single or few dipole solutions can be generated from the magnetic field data 
and used as prior estimates of the current distribution. It is then possible to identify the appropriateness of 
this particular dipole-model prior estimate by computing from Equation 35 the change in the expectation 
associated with including the measurement information without constraining the final solution to a dipolar 
form. Because the statistic provides spatial information it indicates directly those areas where the dipole 
model solution has been modified. Using the variance associated with the a-posteriori current (computed 
from Equation [46|) allows us to plot the number of standard deviations of the change in expectation at each 
point in source space. 

To illustrate the usefulness of this technique, two prior descriptions of the source current for this data 
are postulated - a single moving dipole model and a two moving dipole model. The optimal solutions have 
been found through exhaustive search of the discretised source space by least squares minimisation of the 
fit to the data. The positions of the fitted dipoles for nine time slices are shown in Figure ^. They may be 
compared with the nearest points on the source space to the actual dipole positions. Using each model in 
turn as a prior current and Equations pq and M, the statistical significance of the differences between the 
a-posteriori and the a-priori current distributions are calculated (Figure 0b). 

It can be seen that a single dipole is a good model for the first two time instants and also for the 
ninth. This is not surprising as only one dipole is active at these times. The significance plots for the other 
times suggest systematic discrepancies between the model and the data. The restricted localisation of these 
significant differences points to additional localised sources that have been omitted from the model. The 
next step is the two dipole model (Figure ||c). The two-dipole significant difference diagram suggests that 
this model is adequate for all but four time instants. Comparison with the activation curves identifies these 
as being the times when all three dipoles are active. 

The second illustration uses data from an evoked response study of face-processing using the same exper- 
imental system p3|. Human subjects were presented briefly with photographs of human faces and control 
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Figure 3: a) An isometric projection of the experimental geometry. The dots on the source space are 
the feet of the perpendiculars from the three dipole sources used in this simulation. The activation curves 
(dipole moment vs. time in milliseconds) for these three dipoles are shown in b). The longer period curve 
corresponds to the dominant (central) source. The other curve corresponds to the synchronous and equally 
activated lateral sources. The highlighted detector sites are those for which the data is shown in c) channel 
1 d) channel 2 and e) channel 3. Note that it is only in the case of Channel 3 that the signal follows clearly 
the activation curve for the closest dipole. 



(a) 



(b) 



(c) 



X X 


». 


V 


\' 


' < 


*>v 


h 


'V 


»«. 



• 


• 


f f 

• 


1 1 

• 


!> 


• 


• 


1 

• 


• 



• 
• 


• 
• 


• 


1 

• 


• • 


• 

• 


• 

• 


• 

• 


• 
• 



Figure 4: (a) The maxinium likelihood current density with a zero a priori current distribution. Each frame 
is a two dimensional representation of the shell source space pictured in Figure 2. The 9 frames are snapshots 
at time 0ms, 100ms, 200ms, etc. (see Figure 3). They should be read left to right and top to bottom. The 
optimal regularization parameter was determined by the L-curve method |g^ to be 0.283 x trace{P)/N. The 
crosses on the first frame of this sequence show the positions of the three source dipoles in this representation 
of the source space, (b) The significant differences between a single moving dipole model and the computed 
data. The black dots denote the projected position of the fitted dipole. (c) As (b) but with a two moving 
dipole model. 



objects (e.g. animals) and their neural responses were recorded as a function of time after the stimulus. It is 
known that the early response to face images involves widespread activity in the posterior brain but there is 
limited evidence for the precise distribution and timecourse of the neuronal sources. One suggestion is that 
there are three major areas of activity; in occipital cortex and both right and left ventral occipito-temporal 
cortex ||2j, g5[. Strong occipital activity (starting about 100 ms after the stimulus) is expected to lead to 
concurrent activity in the two other regions with a stronger response in the right hemisphere p3. The 
hypothesised arrangement is therefore similar in geometry to the simulated measurement already discussed. 

Figure || shows the same set of outputs that were presented for the simulated system. In this case, the 
fitted dipoles may be thought of as representing a discrete limited region of source current. Figure oa suggests 
early central activity (frame 3) followed by less prominent localised activity on the right (frame 5). These 
source regions are reflected in the 1-dipole solutions (Figure pb). However, comparison with Figure &) shows 
that the accuracy of the single dipole fit is less than for the simulated data even though the noise levels are 
comparable. This may suggest that there are other active sources present. There is no evidence that these 
are recovered by the 2-dipole model (Figure |t;) as there is little improvement in the significant difference 
maps generated using a two dipole model as the prior (see for example the strong similarity between frame 
4 in Figures Qd and gc). It would be reasonable to infer that the additional sources are diffuse. 

The third illustration relates to the main thrust of the face processing study by Swithenby et al P3| , 
which was to identify statistically significant evidence for differences between the responses to faces and the 
other complex visual stimuli. In the initial analysis the strength of evoked activity within a certain region 
and latency span was parameterised in terms of the signal power integrated over a group of channels and 
a specified latency span. These calculations revealed that the brain activity in the right occipito-temporal 
region following face presentation is significantly different (p=0.05) from activity following non face images 
during the latency span 110 to 170 ms. No other consistent and statistically significant differences were 
found. This data-space analysis, though useful, was complicated by the need to survey the large number of 
possible choices of channel group and latency range. 

The Bayesian framework developed here provides an alternative direct means of directly comparing 
responses to two stimuli. By using one data set as the prior and comparing it with the other data set it is 
possible to identify those regions in source space where there is a statistically significant difference between 
the two source structures within the same source model. We have carried out this calculation for the face 
and control stimuli as a function of position and latency for a simple two-dimensional source space consisting 
of a part spherical shell whose radius is similar to that of the cortical surface (Figure o) . 
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Figure 5: (a) The maximum likelihood current density for the face response with a zero a priori current 
distribution. The 9 frames represent equal steps in time from 70 ms to 241 ms after the stimulus, (b) The 
significant differences between a single moving dipole model and the computed data. The black dots denote 
the projected positions of a single dipole. (c) As (b) but with a two moving dipole model. 



F 




¥ 


4J 







Figure 6: The maximum likelihood current density for the face response with a a priori current distribution 
derived from the control experiment with pictures of animals. The 9 frames represent equal steps in time 
from 70 ms to 241 ms after the stimulus. The regularisation parameter ^ was chosen by using the L-curve 
method [E2|. The optimal value was 1.4 x trace{P)/N. 
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The evidence for statistically significant differences in the right occipito-temporal region at about 155 ms 
is clear. However there is no evidence for differences in earlier latencies, in particular with respect to the 
early source shown in Figure Ha. 

5 Discussion 

The illustrations provided above offer ways of exploiting the new Bayesian results that have been derived. 
Dipolc predictions have been systematically examined using a measure that goes beyond a simple scalar 
goodness of fit measure (i.e. percentage of variance accounted for) to a statistically valid map of fitting 
significance. This addresses a long standing issue, the unreliability of the goodness of fit measure as a 
reliable test of the appropriateness of a given model in explaining a given set of data ||22l. The spatially 
discriminated Bayesian approach gives a test which is more reliable when there is fundamental concern about 
the appropriateness of a model. The second example shows how these ideas may be applied to real data to 
assess the complexity of the dipole model that a given data set can sustain. In similar fashion, the third 
example illustrates how the Bayesian framework allows the direct comparison of two data sets in order to 
identify quantitatively the regions of statistically significant source activities. 

These ideas may be extended further within MEG (and EEG) data analysis. An obvious extension is 
to perform dipole analysis as a precursor to a distributed Bayesian analysis. This may provide a way of 
not only refining information about the depth of source activity but also of assessing the reliability of depth 
estimates. Other possibilities include exploration of the dynamics. For example, times at which there are 
statistically significant change in the data could be identified by using a source distribution estimated from 
each time as the prior for an analysis of the data from the next sampling instant. 

In summary, the analysis presented here comprises a Bayesian estimator of a source current distribution 
in the biomagnetic inverse problem. This is generalisable to other systems and may be used as the engine 
for tests of significant difference in source models and data. 
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